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2 ■ ABSTRACT 

Among the many new irregular sateUites that have been discovered in the 
^ \ last five years, at least six are in the so-called Kozai resonance. Due to solar 

perturbations, the argument of pericenter of a satellite usually precesses from 
to 360°. However, at inclinations higher than ~ 39.3° and lower than ~ 140.7° 
a new kind of behavior occurs for which the argument of pericenter oscillates 
around ±90°. In this work we will concentrate on the orbital history of the 
saturnian satellite S/2000 S5 Kiviuq, one of the satellites currently known to be 
in such resonance 

Kiviuq's orbit is very close to the separatrix of the Kozai resonance. Due to 
perturbations from the other jovian planets, it is expected that orbits near the 
Kozai separatrix may show significant chaotic behavior. This is important be- 
cause chaotic diffusion may transfer orbits from libration to circulation, and vice 
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versa. To identify chaotic orbits we used two well-known methods: the Frequency 
Analysis Method (Laskar 1990) and Maximum Lyapunov Exponents (Benettin 
et al. 1980). Our results show that the Kozai resonance is crossed by a web of 
secondary resonances, whose arguments involve combinations of the argument of 
pericenter, the argument of the Great Inequality (2 A jupUer — 5Xsaturn), longitude 
of the node fl, and other terms related to the secular frequencies g^, qq, and 
Sq. Many test orbits whose precession period is close to the period of the Great 
Inequality (883 yrs), or some of its harmonics, are trapped by these secondary 
resonances, and show significant chaotic behavior. 

Because the Great-Inequahty's period is connected to the semimajor axes of 
Jupiter and Saturn, and because the positions of the jovian planets have likely 
changed since their formation (Malhotra 1995), the phase-space location of these 
secondary resonances should have been different in the past. By simulating the 
effect of planetary migration, we show that a mechanism of sweeping secondary 
resonances, similar to the one studied by Ferraz-Mello et al. (1998) for the aster- 
oids in the 2:1 mean motion resonance with Jupiter, could significantly deplete 
a primordial population of Kozai resonators and push several circulators near 
the Kozai separatrix. This mechanism is not limited to Kiviuq's region, and 
could have worked to destabilize any initial population of satellites in the Kozai 
resonance around Saturn and Jupiter. 

Subject headings: celestial mechanics, satellites, resonances 

1. Introduction 

Garruba et al. (2002) studied an analytical model of the Kozai resonance (Kozai 1962) 
based on a secular development of the disturbing function, expanded in scries oi a/a' (where 
the ' stands hereafter for quantities related to the perturber), truncated at second order 
in a/a', and averaged over the mean anomalies of both perturber and perturbee (see also 
Innanen et al. 1997). Two kinds of behavior are possible in this simplified model: for 
inclinations of less than 39.23° (or higher than 140.77° for retrograde satellites), the argument 
of pericenter freely circulates from 0° to 360°, while at intermediate inchnations a new class 
of solutions, in which the argument of pericenter librates around ±90°, are possible. Such 
behavior is called the Kozai resonance. The value of the critical inchnation of 39.23° is, in 
fact, an artifact of the averaging; the real boundary between the region in which libration 
is possible or not is actually a more complicated function of the satellite's semimajor axis 
and eccentricity (Cuk et al. 2003, Cuk and Burns 2004). Gurrently, six irregular satellites 
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of jovian planets are known numerically to be in the Kozai secular resonance: two around 

Saturn (S/2000 S5 Kiviuq, and S/2000 S6 Ijiraq), two around Jupiter (S/2001 JIO Euporie, 
and S/2003 J20), and two around Neptune (S/2002 N2 and S/2002 N4, Holman et al. 2004). 
A third neptunian satellite, S/2003 Nl is on a chaotic orbit that may also be in the Kozai 
resonance (Holman et al. 2004). 

A question that Carruba et al. left unanswered regarded the possible presence of chaos 
and chaotic diffusion at the separatrix between the regions of circulation and libration. 
Since in that work the authors considered a three-body problem, but averaged over the mean 
anomalies of both pcrturbcr and pcrturbcc, our analytical model was an integrablc onc-dcgrce 
of freedom model. In the real system, however, perturbations from other jovian planets and 
short-period terms may produce chaos near the Kozai separatrix. Chaotic evolution might 
lead to the escape of objects whose orbits were originally inside the Kozai resonance. Here we 
test the hypothesis that the chaotic layer is due to the fact that the Kozai resonance overlaps 
with secondary resonances near the separatrix. It is well known that resonance overlap is one 
of the major sources of chaos in dynamical systems (Chirikov 1979). Moreover, the positions 
of the planets may have changed after their formation, due to gravitational scattering of 
planetesimals (Malhotra 1995). As a consequence of the different planetary positions, the 
shape of the chaotic layer and the locations of the secondary resonances inside and outside the 
Kozai resonance should have been different in the past. This might have had consequences 
on the stability of any primordial population of Kozai resonators. 

This paper will concentrate on identifying the chaotic layer at the transition between 
circulation and libration for Kozai resonators, in particular for the saturnian satellite Kiviuq, 
whose orbit is very close to the separatrix of the Kozai resonance. Our goal is to understand 
the origins of chaos and the effects of planetary migration. 

Our investigation is divided in the following way: section 2 describes the major features 
of our simplified model of the Kozai resonance and identifies the transition region between 
circulation and libration. Section 3 recalls two methods to identify chaotic behavior in 
dynamical systems, the Frequency Analysis Method and the maximum Lyapunov exponents. 
Section 4 investigates possible causes of the chaotic layer. Section 5 shows how the Kozai 
resonance is crossed by a series of secondary resonances of different strengths. In the final 
section we demonstrate how the effect of planetary migration, combined with the effect of 
the secondary resonances, might have operated to depopulate a possible primordial satellite 
group of Kozai resonators. 
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2. Identifying the Transition between Circulation and Libration: Analytical 

and Numerical Tools 

Carruba et al. (2002) described a simple analytical model of the Kozai resonance. This 
model is based on a development of the disturbing function in a series of r/r', where r and r' 
are the radial planetocentric distances of the satellite and the Sun, truncated at the second 
order. The perturbing function R is: 




where Go is the gravitational constant, m' is the solar mass, s is the angle between directions 
to the sateUite and the Sun as seen from the planet, and P2 is the second-order Legendre 
polynomial. The perturbing function, when averaged over the mean anomalies of both Sun 
and satellite, becomes: 



R = [2 + Se^ - (3 + 12e2 - ISe^cos^ u) sin^ jl , (2) 

where a, e, and / are the semimajor axis, eccentricity, and inclination of the satellite, respec- 
tively, b' = d\J\ — e'^ is the Sun's semiminor axis. The perturbing function, when expressed 
in terms of Delaunay elements (/, gi, /i, L, G, H\ see Murray and Dermott 1999, Danby 1988), 
is: 



R 
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(3) 



Since both / and h are cyclic, the conjugate quantities L and H are constants of the motion, 
as is the quantity: 



e= {1 - e') cos' I = HyL'. (4) 

With L and H fixed, the perturbing function (3) represents a one-degree of freedom model for 
the Kozai resonance (Innanen et al. 1997). The following equations for the orbital elements 
hold: 



dl 



15 

levT^e^ 
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— eV 1 — sin 2a; sin^ /, (6) 
8 

, ^ [2(1 - e^) + 5sin2 uj(e^ - sin^ /)] , (7) 

where r = Goui't/nb'^ and t is time. The time has been rescaled in this way to ehminate the 
equations' dependency on the sateUite's semimajor axis. Kinoshita and Nakai (1999) found 
the solution of (5-7) in terms of elUptic functions. Fig. 1 shows the results of numerical 
integrations of equations (5-7) for 6 = 0.70 and 6 — 0.25. For O less than 0.6, librating 
solutions are possible. 

This model has several limitations. First, short-period terms have been eliminated by 

averaging. One consequence of this averaging is that the limit between the zones in which 
libration is possible or not is / = 39.23° (or 140.77° for retrograde orbits), regardless of 
the satellite's semimajor axis (both inclinations correspond to = 0.6 for e = 0). When 
short-period terms are included, the boundary shifts slightly (Cuk and Burns 2004). A case 
where this is observed is the orbit of the jovian satellite Euporie that is a Kozai resonator 
but has e > 0.6. 

A second limitation of the averaged model is that, since it is an integrable, one-degree- 
of-freedom model, no chaos can appear at the separatrix between circulation and libration. 
Finally, since mutual perturbations among jovian planets are ignored, the orbital elements 
of the planet about which the satellites circle are constant. In the real solar system, the 
eccentricities of Jupiter and Saturn vary with several prominent periods associated with 
the secular frequencies and (timescales of up to l/g's ~ 300,000 yr), Great-Inequality 
terms 2Aj — SA^, where the suffices J and S hereafter stand for Jupiter and Saturn, (with a 
timescale of 883 yr), and short-period frequencies, connected with the orbital periods of the 
planets (12.5 and 29.7 yrs for Jupiter and Saturn, respectively). The shape of the separatrix 
may change in time according to the variations of the planet's eccentricity and semimajor 
axis (remember the presence of the factor b'^ in r) which might introduce chaos along the 
real system's separatrix. 

To understand the behavior of the real system, we concentrate on the case of Kiviuq, 
a saturnian satellite in Kozai resonance whose orbit is very close to the separatrix between 
libration and circulation in u. We start by numerically searching for the transition between 
circulation and libration. To do so, we used the following procedure: we numerically inte- 
grated Kiviuq's orbit over 100,000 yrs and found the maximum eccentricity, because then the 
satellite is closest to the separatrix. We computed © according to the inchnation of Kiviuq 
at that time, and generated a grid of initial conditions for test particles; this contains 19 
values in a;, from 40° to 140°, separated by 5°. All particles were given the average value of 
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Kiviuq's semimajor axis during the integration (a = 0.0755 AU), and a grid in inclination 
with 21 values, starting at 38.6° with a step of 0.6°. The values of the eccentricities were 
computed from B = 0.4464. Initial Q and M (mean anomaly) are equal to those of Kiviuq 
at its maximum eccentricity. This grid forms our "low- resolution survey" . 

These initial orbits were integrated for 2 Myr under the influence of the four jovian 

planets (Outer Solar System, OSS hereafter) with a modifled version of SWIFT- WHM, the 
integrator that uses the Wisdom and Holman (1991) mapping in the SWIFT package (Levi- 
son and Duncan 2000). We modified the integrator so that the planets' orbits arc integrated 
in the heliocentric system and the particles' orbits are integrated in the planctoccntric ref- 
erence frame. This setup minimizes the integration error (Nesvorny et al. 2003). 

To more easily identify secular oscillations in Saturn's inclination, we refer all the po- 
sitions to Saturn's initial orbital plane. Fig. 2 shows the results of our simulations for our 
low-resolution survey. The actual separation between circulating and librating behavior does 
not follow the solution of the secular model when the full perturbations from jovian planets 
are considered. Moreover, some orbits alternate between circulation and libration (and, in a 
few cases, spend some time in the other libration island around uj — 270°). 

The question that remains to be answered is whether there is a chaotic layer near the 
separatrix, and if chaotic diffusion might be effective in depleting a primordial population 
of Kozai resonators. To address this question, and identify chaos, we will use the Frequency 
Analysis Method (FAM, hereafter; Laskar, 1990), and will compute the Maximum Lyapunov 
Exponents (MLE hereafter; Benettin et al., 1980) which we now describe. 

3. Numerical Tools to Study Chaotic Orbits 

The Frequency Analysis Method (FAM) is essentially an analysis of the frequency power 
spectrum of appropriate combinations of numerically determined orbital elements. In our 
case, for example, we are interested in the frequency related to the precession of the argument 
of the pericenter uj. Fig. 3 shows the periods in uj for the particles in our low- resolution survey. 
The periods peak sharply in the transition region from circulation to libration in u, reaching 
a maximum of 1800 yrs. The minimum libration period is 480 yrs. The average period of 
libration in the Kozai resonance in our survey is 640 yrs. This range in periods corresponds 
to a frequency range of ~ 700-3200 "/yr. 

We took the Fourier transform of e ■ exp uo — e cos uj + lc sin uj (where u — -\/— T) over 
two time intervals (0-1 Myr and 1-2 Myr) and computed the quantity: 



- 7- 



(8) 

where is the frequency over the first interval and p the corresponding frequency on the 
second interval. If the motion is regular, f{t) will be nearly constant, with small variations 
produced by errors in the determination of the frequencies (to be discussed later in this 
section). If the motion is chaotic, variations in a will be related to changes in a;'s precession 
period. 

While this is the main idea behind the FAM, a few complications occur. First, the 
precision of a Fast Fourier Transform (FFT hereafter) is limited by its coarse resolution 
[A/ = jj^, where At is the (constant) time separation of cj-sampling and N is the number 
of data points used, which must be a power of 2 (Press et al. 1996)] in recovering the 
frequencies having the largest amplitude. To overcome this problem, we used the Frequency 
Modified Fourier Transform method (FMFT), with quadratic corrections (Sidhchovsky and 
Nesvorny 1997), which allows us to retrieve the frequencies with the largest amphtude to 
a better precision. To estimate the error associated with this method, we computed the 
frequencies with and without quadratic corrections, and computed the average value of the 
differences in a. This procedure gives an upper limit on the error. In our case, the average 
error of FMFT was a = —4.5, while a simple FFT had a resolution of o" = —3 (we used a 
time step Ar=30 yr, and 32,768 data points, so that each interval for the determination of 
the frequencies was ~ 1 Myr). 

Another problem to take into account is aliasing. The FFT can recover frequencies 
up to the Nyquist frequency (/tv = ^^^^ corresponds to two samplings per period). 

Frequencies higher than the Nyquist frequency are not lost, but are mapped into a smaller 
apparent value of /, given by fapp = Jn — {freai ~ Jn)- This phenomenon is called aliasing 
and can be easily understood by considering a sinusoidal wave, sampled at isolated points: 
at least two points per period are needed to obtain a good estimate of the period . If 
the sampling is longer than one period, than the retrieved period is larger, and its "ghost" 
signature may appear in the frequency range of interest. 

To avoid this problem, we used an "on-line" low-pass filter, following the procedure of 
Quinn et al. (1991) (based on the work of Carpino et al. 1987). Our filter suppresses to a 
level of 10~^ all Fourier terms with periods smaller than 66.7 yr (stop-band), and attenuates 
to a level of 10~^ all Fourier terms with periods between 66.7 and 200 yr. Terms with periods 
larger than 200 yrs are in the passband (see added information on website). By following this 
procedure, we have eliminated all frequencies connected with the orbital periods of Saturn 
and the satellites. 
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Finally, the shape of the frequency peaks can be modeled as Gaussian, but our mea- 
surements are made at isolated points. This introduces what in signal processing is called 
"impulse noise", which can be attenuated by the use of a median filter (Press et al. 1996, 
Pitas 2000). More details about the use of the median filter and of the test we devised 
to prove the efficacy of our filter can be found in the Appendix. 

To find an alternate way to distinguish between regular and chaotic behavior, we also 
computed Maximum Lyapunov Exponents (MLE) for our set of initial conditions. A detailed 
explanation of the theory of Lyapunov exponents goes beyond the scope of this paper; instead 
we refer the reader to Lyapunov (1907) and Benettin eA, al. (1980). MLE is a measure of 
exponential stretching of nearby orbits. The Lyapunov exponents are equal to zero for regular 
orbits (they tend to zero in finite-time calculations), while they assume positive values for 
chaotic orbits. The inverse of a Lyapunov exponent is the Lyapunov time T^. Smaller values 
of Tl indicate enhanced local stochasticity. 

To estimate MLE for orbits we used a modified version of SWIFT-LYAP2, a code that 
integrates the difference equation (MorbideUi 2002) in the SWIFT package (Levison and 
Duncan 1994). The code was modified to reduce the integration error by integrating the 
planets in the heliocentric frame and the satellites in the planetocentric frame. For each 
of the test particles, we integrated the difference equation with an initial difference vector 
of modulus (i(0) = a/6 ■ 10~^ and determined the modulus d(t) of the displacement vector 
between the two vectors at time t. We constructed a series [t. In [(i(t)/(i(0)] and performed 
a linear least-squares fit on this scries. Since d{t) ~ (i(0) exp {Lt), where L is the Lyapunov 
exponent, the slope of In [d{t) / (i(0)] versus time is equal to the Maximum Lyapunov exponent. 
We computed the Lyapunov exponents for all orbits in our low-resolution survey. We will 
show these results in the following section. 

We should point out that FAM and MLE, while both useful tools to investigate chaotic 
behavior, measure different things. FAM measures macroscopic changes in frequencies (i.e., 
the speed at which chaotic motion changes the frequencies), while measures the rate 
of exponential stretching of nearby trajectories. These two techniques are complementary: 
FAM is more practical in finding regions that are macroscopically unstable, and Tl can reveal 
the presence of different chaos-generating mechanisms. We will discuss in more detail these 
aspects of the two methods in the following section. 
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4. Origins of Chaos in the Kozai Resonance 

To investigate the presence of chaos, we apphed the FAM to our low-resolution survey, 
integrated under the influence of the four jovian planets. Fig. 4 plots logo" for our low- 
resolution survey when the time interval between the first and the second frequencies was 
about 1 Myr. Fig. 4 shows that high values of a are found near the actual separatrix. The 
chaotic layer is asymmetric between values of a; > 90° versus uj < 90°, and Kiviuq seems to 
be on the border of the chaotic layer. 

Having shown that a chaotic layer exists, the next logical step is to ascertain its origin(s). 
Let us recall the analytical model presented in Section 2. In that model, the shape of the 
separatrix depends on Saturn's semiminor axis 6', which is a function of the semimajor axis 
and the eccentricity. Variations in Saturn's semimajor axis and eccentricity may therefore 
generate complex orbital evolution near the separatrix. But on what time-scales? And how 
can we prove this hypothesis? 

To study whether a variable h' may generate chaos, we performed a numerical simulation 
with the Sun (Saturn) on an orbit of fixed eccentricity (equal to the average value of the 
saturnian eccentricity) and the set of test particles used for our low-resolution survey. The 
difference compared to the integrable model is the presence of short-period terms. We 
performed the FAM on test orbits in our low-resolution survey and found that all orbits are 
regular in this case. In contrast, when only Saturn and the Sun were present, and Saturn 
had a constant eccentricity equal to the mean value computed during the 46,000 yr period 
of oscillation (= l/^fe), no chaos was observed. This seems to confirm our intuition that 
variations in Saturn's eccentricity and semimajor axis may be responsible for the chaotic 
layer. 

To test this hypothesis, we performed a simulation in which both Saturn and Jupiter 
are present. In this case the layer of chaos appears (Fig. 5). The similarity between Figs. 4 
and 5 suggest that the effects of Uranus and Neptune are unimportant. To demonstrate 
this, we performed two simulations with our low-resolution survey, one with the complete 
outer solar system (OSS), and the other with just Jupiter and Saturn. To show that the 
results are nearly equal we applied the test (see Appendix) to the filtered results of both 
simulations. We obtained a value of of 22.4 (out of a maximum number of 399) indicating 
a very high probability that the two distributions are the same. While perturbations from 
Uranus and Neptune might play a role in modifying the chaotic layer, such perturbations 
are clearly less important than Jovian perturbations. So, in the following, we only consider 
jovian perturbations. In doing so we of course modify the values of planetary frequencies 
and mean motions, and therefore the positions of any secondary resonances. But, since in 
this phase of our work we are interested in understanding on what timescales the chaotic 
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behavior is driven, we believe this is an acceptable approximation. 

There are two prominent timescales for oscillations of Saturn's e and a: the secular 
oscillations (period of order l/grg ~ 46,000 yrs) and the oscillations related to the Great 
Inequahty (883 yrs). To study which of these perturbations is more relevant to the ori- 
gin of chaos, we describe the OSS' orbital evolution, in terms of Fourier series (Bretagnon 
and Francou 1988). In this model, the equinoctial orbital elements {a,X,K — e cos zu,H = 
e sin tu, Q = sin | cos Q,P = sin | sin Q) are developed in series of cosines and sines of combi- 
nations of twelve angles: four mean longitudes of the giant planets, four proper longitudes of 
pericenter, and four proper longitudes of the nodes. These proper angles are linear functions 
of time. This model assures a fractional precision better than 5 ■ 10~^ for K, H, Q, and P 
over 10 Myr. 

To investigate what minimal model can create a chaotic layer similar to the one observed 
in the full problem, we used different solutions for the orbits of both Jupiter and Saturn. We 
modified our version of SWIFT- WHM so that both orbital evolutions of Jupiter and Saturn 
are computed "on-line" according to the Bretagnon secular solution. The test particles are 
then subjected to perturbations of planets that evolve on such orbits. We considered four 
models for Saturn's orbital evolution: (1) SEC, in which we limit the expansion to terms 
containing the g^.ga, and frequencies for K,H,Q, and P, assume a constant semimajor 
axis, and take 12 terms for A's development. (2) A21SEC, in which we have included 12 terms 
associated with the 1:2 commensurability between Jupiter's and Saturn's mean motions, and 
terms associated with the Great Inequality for the development in a, plus the same terms as 
in the SEG solution. (3) GISEG, in which there were 8 terms for K and H, and 6 for P and 
Q associated with the Great Inequality, and a constant semimajor axis; and (4) A21GISEG, 
which includes the 2:1 terms in a, the Great-Inequality terms and the secular variations. 
Table 1 gives the number of terms taken in a. A, {K, H), (Q, P) for each solution. 

Fig. 6 shows the results of the integration with these models. All simulations used the 
same solution for Jupiter (10 terms in the development of a, 2 terms for K and H, 2 for Q 
and P, and 12 for A). These figures may be compared with Fig. 5, where the traditional 
version of SWIFT- WHM was used. 

Table 2 presents the results of a test in which we compared simulations with "Bre- 
tagnon" models with those of the full integration of Jupiter and Saturn. According to these 
results, terms associated with the Great Inequality and the 2:1 terms in the development of 
Saturn's semimajor axis are the most important for the presence of the chaotic layer. The 
major sources for chaos are the variations in Saturn's eccentricity connected with the Great 
Inequality. 
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To further investigate this hypothesis, we computed the maximum Lyapunov exponent 
for the particles of our low-resolution survey, in a model containing Jupiter and Saturn. 
Fig. 7 shows the resulting values of (left) when the full problem is considered and (right) 
when the planets evolve according to the SEC solution. A remarkable asset of Fig. 7a is the 
region of high values of Tl (i.e., a smaller degree of local stochasticity) near the vertex of the 
separatrix, which, according to the FAM, is a region of high dispersion of the frequencies. 
In this region ^ 1, 700 yrs, i.e., twice the period of the Great Inequality; thus our results 
suggest that the feature is connected to a secondary resonance with the Great Inequality. 
Note, however, that all orbits seem to have quite small values of Tl. Orbits inside the 
libration island, which, according to FAM, are macroscopically regular, have Tl < 2, 000 
yrs. This appears to be a general characteristic of high-inclination irregular satellite orbits, 
as all seem to have Tl < 10, 000 yrs (M. Holman, private communication, 2003). 

The existence of a region characterized by relatively long Lyapunov times, on one hand, 
but high frequency dispersion, on the other, seems counter-intuitive at first. However, recall 
that the rate of frequency dispersion is a macroscopic characteristic of the orbit that is 
related to the global structure of the phase space region. In contrast, the Lyapunov time is 
a local property whose value is dictated by the mechanism that is responsible for the chaotic 
motion. 

The existence of regions of large Tl but rapid diffusion is not uncommon in the asteroid 
belt. For example, in the 7/3 resonance (sec Tsiganis et ai, 2003), the smaller values of Tl 
arc observed at the borders of the resonance where they arc related to the pulsation of the 
separatrix due to secular precession. However, the macroscopically most unstable region is 
located at smaller libration amplitudes, is characterized by larger values of Tl, and is related 
to the colocation of the uq secular resonance. Similar results can be shown for other low-order 
mean-motion resonances (e.g., 2:1 and 3:1), where the most unstable zones are generated by 
coexisting secular resonances. As T^'s value is typically of the order of the forcing period 
(e.g., 1/5^5 305,000 yrs, l/ge ^ 46,000 yrs, Tsiganis and MorbideUi 2004), these regions 
generally have higher values of Tl than the region associated with the pulsating separatrix 
of the mean-motion resonance. 

We are likely observing a similar phenomenon in the case studied here. We believe that 
this high-Tj, region is associated with a secondary resonance, involving the Great-Inequality 
period. This is also supported by the results shown in Fig. 7b (SEC solution), where all orbits 
inside the libration zone become regular (note the different scale compared to Fig. 7a). To 
further support our claim, we selected a particle from this high-T^ region, for which we plot 
the time evolution of the argument 'nj — fl — 5Xs + 2Xj + 3Q.s (Fig. 8). This corresponds to one 
of the possible critical arguments of the 1:1 resonance between the frequency oi cu — zu — ^ 
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and the frequency of the Great Inequahty, 2Aj — 5A5 (3^5 is added in order to fulfill the d' 
Alembert rules of permissible arguments, i.e., the sum of coefficients of individual longitudes 
must be zero and the sum of coefficients of nodal longitudes be even). The behavior of the 
arguments switches erratically between intervals of libration and circulation, which indicates 
a transition through the separatrix of this secondary resonance. 

The following picture arises from our low-rcsohition surveys. Two major sources of chaos 
exist in the system. One is connected to the transition region from circulation to libration 
and, in particular, is driven by variations in Saturn's eccentricity, mainly associated with 
Great-Inequality terms, and variations in Saturn's semimajor axis. The other is related to a 
secondary resonance involving the argument of pericenter uj and the Great-Inequality terms 
for circulating particles. In most cases this resonance is a stronger source of chaos than the 
transition region. Better resolution is needed to understand the behavior and shape of this 
resonance (see Sec. 5). 

One problem that remains to be explained is the asymmetry of the chaotic layer with 
respect to a; = 90° (see Figs. 5-7). According to the secular model, the result should be 
symmetric. Our simulations with the Bretagnon model suggest that variations in Saturn's 
semimajor axis may be partly responsible for the asymmetry. Alternatively, the asymmetry 
in the chaotic layer may be related to our choice of initial conditions. To exclude this 
second possibility, we need to ascertain that our choice of initial conditions {Vt and M) 
for the test particles' orbits is not introducing spurious effects. To our knowledge, three 
resonant configurations may alter the eccentricities of the test particles, and so put their 
orbits closer or farther away from the region of chaos that we are interested in. The three 
resonances are: "nodal" (Q — Jl© ^constant), pericentric {vo — vdq — 0), and the evection 
inequality ^{vo — Xq) ^constant (= 45°, so as to have an intermediate value of the evection 
angle)^. In our simulation, the Sun's orbit starts with: Q, — 23.162°, a; = 311.992°, and 
M — 322.291°. We computed Q, for the test orbits according to the resonant configuration. 
For the pericentric resonance, we used two values of the resonant argument: w — Wq = 
and (Yokoyama et ai, 2003), w — vjq = 180°, X — Xq = 180°. The second configuration 
maximizes the perturbing effect of the pericentric resonance. 

To have a quantitative measure of the asymmetry, we use an "asymmetry coefficient" 
that is the difference of cr's values for uj larger and smaller than 90°. To compute this 
coefficient for orbits in the libration island (60° < uj < 120°), we take two columns symmetric 
about LU = 90° (e.g., the columns for uj = 60° and uj = 120°, with the obvious exclusion of 



^ Since the evection term in our case is not resonant (Nesvorny et al. 2003) we refer to the angle 2{w — Xq) 
as evection inequahty. 
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the column for uj — 90°) and compute the fractional difference for each pair. We repeat 
the process for all pairs of values and compute the average and standard deviation of the 
measurements. The average value is our "asymmetry coefficient" , and the standard deviation 
gives an estimate of the error. Obviously, the lower is the asymmetry coefficient, the more 
symmetric is the distribution of a around u) = 90°. 

Table 3 gives asymmetry coefficients for each simulation. Since the secular resonance 
with resonant argument equal to 0° seems to give the more symmetric configuration, we 
believe its effect might be the dominant one (but large errors prevent our firm conclusion). 

5. A Web of Secondeiry Resonance Crosses the Kozai Resonance 

The low-resolution survey has yielded information regarding the large-scale structure of 
chaos inside the Kozai resonance. However, because most of the chaotic behavior is present 
near the boundary region between circulation and libration, we now concentrate on analyzing 
this region in more detail. For this purpose, wc constructed a new set of initial conditions: 
we used 276 bins in inclination, starting at 37.5° and separated by 0.02°. We used one 
UJ (= 90.0°). The choice of the other orbital elements followed the same criteria used for 
the low-resolution survey. We called this new set of initial conditions a "high-resolution" 
survey. We concentrate on the region around u — 90° to eliminate the problems associated 
with asymmetries in uj suggested by our low-resolution survey. Results of integrations of 
test particles with uo equal to 89.8° or 90.2° are very similar to those with u — 90° and are 
not shown. Our choice of values for the particles' inclinations allow us to sample periods in 
libration and circulation that go from a minimum of 450 yrs to a maximum of 820 yrs for 
librating particles, and from a minimum of 650 yrs to a maximum of 1820 yrs for circulating 
particles, thus covering the whole transition region. 

We integrated the test particles of our high-resolution survey under the influence of the 
four jovian planets for 2 Myr and applied FAM. Fig. 9 shows two regions of high chaoticity. 
The one at x ~ 0.73 is clearly associated with the transition from circulation to libration: 
test particles switch from a;-circulation to a;-libration at this x (dots in the plot). The 
other region of strong chaotic behavior at x ~ 0.68 is associated with the pericentric secular 
resonance {vo — — 0). Fig. 10a shows the time evolution of the resonant argument for an 
orbit in this resonance (xq = 0.6848). 

More interesting for our purposes is the region at x ~ 0.696. A plot of the resonant 
argument ro — Q — 5A5-|-2Aj-|-3Q5 shows that this feature is associated with the 1:1 resonance 
between uj and the Great Inequality (other resonances involving wg and wj^ instead of Q5, 
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are also present but are weaker. This property is generally shared by other resonances 
involving harmonics of the Great-Inequality period, like the 4:3, 3:2, etc.). The resonance we 
discovered is only present for circulating particles. A question that might arise is if a similar 
resonance could also be found for librating test particles. To answer this question we consider 
the periods of libration and circulation and estimate them by the frequency associated with 
a;'s precession found via FAM. Fig. 11 shows such periods as a function of x. We also report 
the regions having high values of cr, possibly associated with secondary resonances (vertical 
lines), and the commensurabilities between periods in uj and Great- Inequality periods (for 
example, the 3:4 commensurability means that, for three periods in uj there are four periods 
in 2Aj — 5A5). 

We are conscious that the commensurabilities between the Great Inequality and uj are 
not real resonances, since they do not respect the d' Alembert rules. Nevertheless, we think 
this is a useful plot, since gives a first-order estimate of the position of the real resonances. 
To compute the expected positions of the actual resonances we used the following procedure: 
for a resonance of argument vo — Vt — b\s + 2Aj -|- 3O5 (we will use a; for — O hereafter), 
the frequency of a;'s precession for which there is resonant behavior is given by: 

/a; = /G/-3^, (9) 

where fci is the frequency associated with the Great Inequality term (2Aj — 5A5, equal to 
1467.72 "/yr for the current position of the planets), /^g is the frequency associated with 
(= -26.345" /yr, Bretagnon and Prancou 1988, the same source is used for the values of 
and go), and is the precession frequency. Analogous equations can be used for other 
resonances. From a;'s frequency it is straightforward to determine the period of precession, 
and from that period we can determine the value of x. Fig. 11 seems to be very instructive. 
For example, since the period of the 1:1 resonance is 840 yrs and since the first test particle 
in librating regime has a period of 820 yrs, at the end of the transition region, there is 
currently no equivalent of the Great-Inequality secondary resonance for librating particles. 
Also, apart from the two regions with a associated with the pericentric secular resonance and 
the Great-Inequahty resonance, another interesting feature of high chaotic behavior, which 
seems to be associated with a strong secondary resonance, appears at x — 0.73. Fig. 11 
shows how in this region the values of periods in uj are constant, which is an indication of 
librating behavior in a secondary resonance. Other regions of high a values can be found 
at X = 0.722 and, in the librating region, at a; = 0.741 (since this region is so close to the 
separatrix, the chaotic behavior here observed might be related to the pulsating behavior of 
the separatrix itself, when perturbations from other jovian planets are considered). 

We beheve that the features of high a x — 0.722 and x — 0.73 are associated with 
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two other secondary resonances, whose resonant arguments are given by A{zu — Q) — 3(5A5 — 
2Xj)+8Qs+'^s and 3(ro — Q — 5A5 + 2Aj + 3Q5) + (Q — respectively. Fig. 10c and d show 
the behavior of the argument for orbits near these two resonances. Other weaker resonances 
connected with other commensurabihties between u and 2Aj — 5A5 are also observed. The 
resonance of argument 3{zu — fl) — 2{5\s — 2Aj) + 5Qs + '^s is expected to be stronger 
than the one associated with the 3:4 commensurability, but its location is so close to the 
separatrix that its effect is difficult to discern. Table 4 hsts the resonances for which the 
librating behavior of the resonant argument is observed, and a few candidates that could 
explain features of weaker chaos. 

So far wc have only discussed the case of circulating test particles. Features of weaker 
chaos, however, exist also in the librating region. Unfortunately, in this case plotting the 
resonant argument is not so easy, since, by the definition of libration, cu oscillates around 
±90° instead of covering all values from to 360°. A possible way to overcome this problem 
would be to make a change of variables, so as to put the origin of the system of coordinates 
for ecoscj and e sin 00 at the libration center. The new angle 00' would then rotate from 
to 360°, and it would then be possible to check for the behavior of the resonant arguments 
by combining uj' with Gl and other terms. This procedure is rather cumbersome, especially 
since the libration point itself is not a fixed point, but, due to perturbations from the other 
jovian planets, oscillates with timescales associated with the Great Inequality, g^, qq, se, etc. 
Fortunately, there is an alternative method to determine whether a resonance has in its 
resonant argument terms associated with the Great Inequality. 

It is widely believed that planets have migrated since their formation (Malhotra 1995, 
Gomes 2003, Levison and Morbidelli 2003). In particular, the gravitational scattering of a 
planetesimal disk modified the positions of the jovian planets so that, while Saturn, Uranus 
and Neptune migrated outward in semimajor axis, Jupiter, being the most effective scatterer, 
migrated inward. As the work of Malhotra suggested, the origin of the highly eccentric, 
inclined, and Neptune-resonance-locked orbits of Pluto and the Plutinos might be explained 
in the context of sweeping resonant capture due to the changing position of Neptune's orbit. 
Our work has shown that a major source of chaos for test particles in Kozai resonance around 
Saturn is due to a secondary resonance between the argument of pericenter lv, the argument 
of the Great- InequaUty (i.e., 2Aj — 5A5), and terms connected with the planetary frequencies 
g5,ge, and sq. It is well known that by altering the positions of Jupiter and Saturn, the 
period of the Great Inequality changes. Fig. 12 shows how, by varying the initial value of 
the osculating semimajor axis of Jupiter by a positive amount (and keeping the position of 
Saturn fixed) it is possible to modify the period of the Great Inequality from the current 
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value of 883 years to values of ~ 400 yrs or less. ^ 

The average and minimum values of a;'s libration periods for particles in Kiviuq's region 
are 640 yrs and 480 yrs, respectively. Presently, the Great-Inequality resonance has one 
location, just outside the separatrix, in the region of circulating particles. An interesting 
question would be to see what happened when the Great-Inequality resonance had a shorter 
period. In particular, when the value of the Great Inequality was shorter than ~ 820 yrs 
(the maximum period of libration), a new secondary resonance should appear in the region 
of librating orbits. When the Great-Inequality's period was 640 yrs, several particles in the 
libration island would have been in the region of this new secondary resonance. This can 
be seen in Fig. 13 which shows the results of the integration of our low-resolution survey 
integrated under the influence of Saturn and a Jupiter whose orbit was modifled so that the 
period of the Great Inequahty was a) 640 yrs, and b) 480 yrs. We called these simulations 
"static integrations", as opposed to other simulations to be discussed in the next section 
for which the period of the Great Inequality is not fixed but changes with time. Note how 
the position of the high-chaoticity region, that our simulations show to be connected to the 
1:1 resonance between uj and the Great Inequality, is displaced upward with respect to the 
traditional integration for the first case, and is very close to the center of the libration island 
for the second case. More important, a significant fraction (25%) of orbits that stayed inside 
the libration island in the integration with the present configuration of the planets were on 
switching orbits when the period of the Great Inequality was 640 yrs. 

These simulations show that indeed a 1:1 resonance between u and the GI (and the 
other terms needed to satisfy the d'Alembert rules) was present in the libration region when 
the GI's period was lower. Goming back to our problem of identifying secondary resonances 
for the high-resolution survey, the method applied to our low-resolution survey to change the 
GI's period can give precious insights on the identity of those resonances. Fig. 14 shows how 
the frequency with largest amplitude in the spectra obtained with FAM (just the frequency 
with largest amplitude, not necessarily the frequency associated with a;'s precession) for our 
high-resolution survey changes with x. When a secondary resonance is encountered, values 
of the frequencies spread around rather than following a quasi-linear behavior (the feature 
associated with the pericentric secular resonance at a; = 0.68 is most instructive). 

To test our interpretation of the secondary resonances, we modified Jupiter's position so 



similar change can also be achieved by modifying the initial value of the difference in mean longitude 
of Jupiter and Saturn (Ferraz-Mello et al. 1998). This alternative method has the advantage of better 
preserving the values of the secular frequencies of the planets (i.e., g5,ge,SQ). However, in this work, we are 
interested in the actual effect of planetary migration, which alters the semimajor axes of the planets, and so 
we prefer to use the first method. 
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that the GI period became 810, 883 (current value), and 950 yrs, respectively, and integrated 
the high-resolution survey with the three configurations. Figs. 15a, b, and c illustrate what 
happens to the frequencies for these three values of the GI, for the region of the 1:1 reso- 
nance, the 3:4, and for librating particles, respectively. Fig. 15a shows the behavior in the 
region of the 1:1 resonance. The resonance position (computed using Eq: 9 and the different 
GI periods) shifts toward larger x when the GI period becomes longer. When the period 
equals 810 years, the position of the resonance corresponds to that of the pericentric secular 
resonance (note how its position does not move when the GI period is changed) , and it goes 
to higher values of x when we increase the period. The vertical lines indicate the expected 
position of the resonance when the GI period is modified. The fact that there is an excellent 
agreement between the predicted position of the resonance and the results of our numerical 
simulation seems to further confirm our hypothesis for the source of chaos in this region. 

The middle row of Fig. 15 shows the same plots, but for x — 0.71 — 0.735. This 

region contains two resonances of argument 4:{zu — Jl) — 3(5A5 — 2Aj) -|- 8fls + zus, and 
3{zu — Q — 5Xs + 2Aj -|- 3Q5) -|- (n — Q5): these are given the short-hand notation of 4:3 
and 3(1:1). Once again, the numerical simulations demonstrate that the chaos zones move 
with the resonances, i.e., they agree with the predictions based on our resonant arguments 
(vertical lines). 

Finally, the bottom row of Fig. 15 shows the values of the frequencies of largest amplitude 
for the region of librating satellites. The dashed vertical lines represents the transition region 
in which particles alternate from circulating to librating behavior. The left panel shows 
substantial frequency diffusion near the separatrix, when Pqi = 810 yrs, not observed for 
longer Great- Inequality periods. We believe that this is due to the appearance of the 1:1 
resonance (having argument zu — ft — 5Xs + 2Xj + 3D,s) in the region of librating particles 
(see the vertical line on the plot). This resonance is simply out of the range of periods in cu 
for librating particles when Pqi > 820 yrs. 

Regarding other resonances in the region, we are still limited by the problem of plotting 
the resonant argument. However, our simulations suggest at least two other resonances in 
this region: 2('nj -Q) -3(5Xs -2Xj) +8^3 + '^s and 3('nj -Q) +A{5Xs -2Xj) + 11^3 + '^s- 
The bottom row of Fig. 15 shows how the expected positions of these resonances shift when 
the Great-Inequality period is modified. In each case, the expected position is very close 
to regions of frequency variation. We believe this is strong circumstantial evidence in favor 
of the existence of these resonances. Other resonances containing combinations of only wj 
or zus instead of fls (plus any other term needed to satisfy the second d'Alembert rule) in 
their argument are also possible, but, as observed in other cases, are associated with features 
of weaker chaos (their locations are not shown in Fig. 15). Regarding the bottom row of 
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Fig. 15, we observe that the resonance positions shift toward the separatrix as the GI period 
increases. 

To conclude, we have identified several secondary resonances connected to the Great 
Inequahty for the region near the separatrix of the Kozai resonance. The hkelihood that 
the jovian planets had different past positions means that the locations of these resonances 
might have also been different. This introduces interesting new prospcctivcs for the stability 
of primordial satellites inside the Kozai resonance, that we now investigate. 

6. Effects of Planetary Migration on Any Primordial Population of Satellites 



We have just argued that a web of secondary resonances involving a;, the argument of the 
Great Inequality, and other terms exists in the region of phase space around the separatrix 
between libration and circulation. Considering that the period of the Great Inequality might 
have been different in the past, we ask ourselves if a mechanism of sweeping secondary 
resonances inside the Kozai resonance, not dissimilar from the one that acted in the Kuiper 
Belt, could also have acted inside the Kozai resonance. If that is the case, what are the 
repercussions for the stabihty of sateUites' orbits in the Kozai resonance? 

To address these questions we need to a) have a model of planetary migration and b) 
have an integrator able to simulate the effect of planetary migration for Jupiter and Saturn. 
Following Malhotra (1995), we assumed that the semimajor axis varied as: 



where a/ is the semimajor axis at the current epoch, Aa is the change in semimajor axis 
(equal to -0.2 AU for Jupiter and to 0.8 AU for Saturn) and r (= 2 • 10^) is a characteristic 
timescale for migration. For the integrator to simulate planet migration, we followed this 
recipe. We modified SWIFT- WHM so that an additional drag-force was applied to each 
planet along the direction of orbital velocity (by the symbol v we identify an unit vector 
along the direction of orbital velocity). This produces an acceleration: 
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where v is the acceleration and is the initial position of a planet (Malhotra 1995). We 
neglected the effect of the planetesimals' perturbations on the satellites. Fig. 16 displays 
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several computed evolutions of Saturn's and Jupiter's semimajor axes. The two curves agree 
well. These are to be compared to the results of Beauge et al. (2002) who used SWIFT- 
RMVS3 to simulate the evolution of a disk of 1000 massless planetesimals subjected to the 
gravitational perturbations of the four major outer planets. The evolution of the planet's 
semimajor axes generally follows the exponential law of Eq. 10, but in addition is affected 
by short-time variations due to the impulses of single encounters. 

Having now the tools for our investigation, we must define initial conditions for our 
simulations. To investigate the effect of planetary migration on the stability of orbits in the 
Kozai resonance, we first integrated Kiviuq with the integrator that accounts for planetary 
migration backward in time for 10 Myr. Not surprisingly, Kiviuq remained inside the Kozai 
resonance for the full length of the integration. We then used the last 100,000 yrs of the 
integration to compute the averaged orbital elements of Kiviuq and O's value (= 0.4322). 
Using this information, we generated a low-resolution survey of test particles, exactly as in 
Sec. 2. We then integrated this new set of initial conditions forward in time for 10 Myr. 
During this forward integration, the Great-Inequality's period increased from 90 yrs at the 
beginning of the simulation to the current value of 883 yrs at the end of the integration. 

Fig. 17 shows the particles' fates at left t — {) and right t — b Myr. 14% of the particles 
originally in the libration island were on circulating or switching orbits by t = 5 Myr. At the 
end of the simulation, 15% of the test particles originally in the libration island were no longer 
Kozai resonators. The mechanism of sweeping secondary resonances seems to be effective 
in depopulating the resonance. This mechanism not only affects orbits near the separatrix, 
but also those particles well inside the Kozai resonance. Examples are the particles with ujq 
— 85° and xq — 0.80 and 0.83, which are very close to Kiviuq's orbital region. This shows 
that this region is also affected by the sweeping secondary resonances. 

We explain these results in the following way: when, owing to the planets' migration, 
the secondary resonance's period approaches ~ 400 yrs, a few of the particle inside the 
Kozai resonance (but only a fraction, the process is not 100% efficient) are captured into 
a secondary resonance, and move further out in the libration island until the secondary- 
resonance's period attains ~ 700 yrs. For such libration periods, the orbits are so close to 
the separatrix that escape becomes possible. This mechanism seems to be confirmed by the 
time evolution (Fig. 18) in uj for one of these particles {xq — 0.83, a;o = 85°). 

The fact that most particles escape from the Kozai resonance in the first half of the 
simulation is not surprising. The period of the Great Inequality reaches 600 yrs during the 
first 5 Myr, i.e., the average value of the precession period for orbits in the Kozai resonance 
(see Fig. 3). After 5 Myr, the growing value of the Great-Inequality period sweeps fewer and 
fewer particles inside the Kozai resonance, and the mechanism loses its efficacy. 
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Another interesting consideration is that not only librating particles are pushed toward 
the separatrix, also some circulators in general are captured by one of the secondary res- 
onances and are carried toward the separatrix. The fact that circulators do not generally 
become librators is due to the fact that in order to cross the separatrix they have to reach 
very high values of the precession period in uo (see Fig. 11). Once close to the separatrix, they 
reach a highly chaotic region, and many particles are lost from the secondary resonances. 
This mechanism could explain why so many irregular satellites are currently found in the 
proximity of the Kozai separatrix (Cuk & Burns 2004 for A J) . 

A word of caution should be given about our results. We only used a smooth exponential 
evolution (Malhotra 1995) for the planet's migration. In the real Solar System things were 
most hkely different from this model; e.g., Hahn and Malhotra (2000) introduced short-time 
scale variations in Neptune's outward expansion by adding some random jitter to the torque 
applied to Neptune (cf. the a history by Beauge et al (2002) shown in Fig. 16). This jitter 
was parametrized by the standard deviation ajuter of the planet-migration torque in units 
of the time-averaged torque. For small values of ajuter (< 10), capture efficiencies in the 2:1 
resonances were substantially unaltered, and only when ajiuer reached values larger than 25 
were capture probabilities reduced. 

In our case, we can assume that a similar mechanism should be at work, and that 
any significant amount of short-time variations in the semimajor axes of both Jupiter and 
Saturn might in principle reduce the capture efficiency into the Great-Inequality secondary 
resonance. Since our work is, in many ways, exploratory, we believe the use of the exponential 
model is defensible at this stage of our study. But we acknowledge that a more realistic model 
for the motions of Jupiter's and Saturn's semimajor axes should be used, before realistic 
values of capture probabilities for our mechanism of sweeping secondary resonances could be 
computed. 

In any case, we believe that our mechanism of sweeping secondary resonances might 
be an effective scheme for destabilizing the primordial population of bodies originally in the 
Kozai resonance. The mechanism should not be hmited to Kiviuq's region. The fact that 
the Great-Inequahty's period swept through values starting from 90 yrs at the beginning of 
the simulation with planetary migration to the current 883 yrs, implies that an analogous 
mechanism should also have been at work in the regions of the other Saturnian and Jovian 
satellites presently in the Kozai resonance, for lu libration periods of less than 1000 yrs. 
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7. Summciry and Conclusions 

We have analyzed the orbital behavior of satellites inside and outside the Kozai reso- 
nance for the phase space around the Saturnian satellite Kiviuq. Our goal was to identify 
areas of chaotic behavior and to understand if chaotic diffusion might have played an impor- 
tant role, now or in the past, for orbital stability inside the Kozai resonance. By applying 
the frequency analysis method (FAM) and computing the maximum Lyapunov exponents, 
we identified several secondary resonances. The major source of chaos for orbits in the circu- 
lating region are secondary resonances involving the argument of pericenter uj {= vj — Vt), the 
Great- Inequality's argument (2Aj — 5A5), the longitude of the node (Q), and the planetary 
secular frequencies ((75, g^, and Sg). The region of transition between circulation and libration 
of UJ is strongly chaotic, when perturbations from the other jovian planets are considered. 
According to our secular model, variations in Saturn's semiminor axis (and, therefore in its 
eccentricity and semimajor axis) are what drive the chaotic behavior, especially the trans- 
fer from circulation to libration. While secular variations modulated by the g5,gQ, and Sq 
frequencies are not negligible, once again our results indicate that the dominant effect is 
connected with perturbations having the Great-Inequality's period. 

Also, contrary to what might be expected from the secular model, our results indicate 
that the chaotic layer at the boundary between circulation and libration is not symmetric, but 
shows an asymmetry between values of cu larger and smaller than 90°. This is a consequence 
of other perturbations connected with secular pericentric resonance, evection inequality, and 
variations in Saturn's semimajor axis, which are strongly dependent on the initial conditions 
for particles' w and the Sun's position. Our simulations suggest that the dominant effect 
is connected with the secular resonance, but variations in Saturn's semimajor axis are not 
entirely negligible. 

From our simulations we can introduce the following scenario for chaotic orbits near 
Kiviuq's Kozai resonance. Starting from librating orbits close to the libration center, we 
possibly identify two weak secondary resonances of argument 2(tu — Q) — 3{5\s — 2Aj) + 
8fls + and 3{zu — fl) — 4(5 — 2Aj) + llfls + that introduce chaotic behavior. 
The first strong source of chaos is however connected with the transition region between 
circulation and libration: orbits start to switch back and forth from circulation to libration 
and the whole region is characterized by high values of a. Once in the circulating regime, a 
strong resonance of argument 3(ro — Q — 5A5 -|- 2Aj -|- SQg) + — ^s) appears a,t x — 0.73. 
In addition to that, the dominant source of chaos, excluding the nearby pericentric secular 
resonance, is connected to a secondary resonance with the Great Inequality of argument 
-07 — — 5A5 + 2Aj + 3Qs (these orbits have argument-of-pericenter periods of order 840 
yrs). Orbits in the region of this secondary resonance show some of the highest a values 
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among the test particles we studied. Other secondary resonances connected with the Great 
Inequahty are also observed. 

The fact that the dominant source of chaos for orbits in Kiviuq's region is connected 
with the Great-Inequality secondary resonance opens interesting perspectives. It is believed 
that jovian planets formed in different locations than their current ones, and then migrated 
by gravitationally scattering the primordial population of planetesimals (Malhotra 1995). 
The possibility that planets occupied different semimajor axes in the past has important 
repercussions, since the Great-Inequality's period is connected with the osculating semimajor 
axes of Jupiter and Saturn. In particular, when the two planets were closer together, the 
period of the Great Inequality was shorter. We simulated this effect first by changing the 
position of Jupiter so that the Great- Inequahty's period was 640 yrs and 480 yrs, respectively 
("static integrations"). These simulations showed that, indeed the position of the Great- 
Inequality's secondary resonance shifted higher. A twin resonance for librating particles 
appeared when the Great-Inequality's period equaled the average libration period (640 yrs), 
and the minimum (480 yrs) (this resonance is not visible today because the maximum period 
of libration inside the Kozai resonance is of ^ 820 yrs). Then we fully integrated Jupiter 
and Saturn with an integrator that simulates planet migration according to the exponential 
law of Malhotra (1995). 

Our results show that a mechanism of sweeping resonances, analogous to the one that 
operated in the Kuiper Belt (Malhotra 1995), or inside the 2:1 mean motion resonance with 
Jupiter (Fcrraz-Mello et al. 1998) must have operated as well inside the Kozai resonance. 
Our simulations show that 15% of the satellites originally inside the Kozai resonance became 
circulators or exhibited switching behavior when planets acquired their final orbits. Indeed 
several particles were captured into the secondary resonance with the Great Inequality (or 
other secondary resonances) and, as a consequence, their amplitudes of oscillation in cu 
increased until the test particles reached the boundary between circulation and libration, 
where escape became possible. Several particles on circulating orbits were also pushed toward 
the separatrix by an analogous mechanism of sweeping secondary resonances (but none 
crossed the boundary to become a librator). We believe that a similar mechanism might 
have acted in the past also for the region of phase space associated with the other jovian 
satellites currently inside or near the Kozai resonance. This process could explain why so 
many irregular satellites are currently found near the Kozai separatrix (Cuk and Burns 2004). 

In this work we used several numerical tools. In a few cases, like for example the median 
filter for our FAM results, we introduced a technique, that, while employed in other fields 
(for example in remote-sensing and in the analysis of data from martian probes), was not, 
to our knowledge, used previously in celestial mechanics for reducing FAM results. Still on 
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the technical side, we created new integrators to account for the Bretagnon solutions of the 
giant planets and to simulate planetary migration (Malhotra 1995). We believe these tools 
might be useful for other studies, not only for the irregular satellites of jovian planets. 

Several questions remain that might be the subjects of future research. In this work we 
could not compute reliable capture probabilities for the Great-Inequality resonances when 
planetary migration is considered, since short-period variations in the scmimajor axes of 
Jupiter and Saturn should also be taken into account (Hahn and Malhotra 2000, Beauge et 
al. 2002). It would also be interesting to extend our studies to the neighborhoods of the 
other jovian satellites currently inside the Kozai resonance, and sec if a web of secondary 
resonances (possibly involving the Lesser Inequality, the 2:1 quasi-resonance between Nep- 
tune and Uranus) is also present for the case of those neptunian satellites in Kozai resonance 
(Holman et al. 2004). We hope that this work might have opened interesting new prospects 
for the study of origins and stability of irregular satellites in Kozai resonance. 
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A. Median filter and test 

When measuring frequencies with a FMFT, the shape of the frequency peaks can be 
modeled as Gaussian, but our measurements are made at isolated points. This introduces 
what in signal processing is called "impulse noise" , which can be attenuated by the use of a 
median filter (Press et al. 1996, Pitas 2000). 

In our case we took data on a running 3x3 matrix; the value at the central point is the 
median of the nine elements. Then we moved the matrix by one clement and repeated the 
procedure. To compute values of a for points along the boundary of our figure, we added 
two external rows and columns of points all with values of a obtained from the average 
of all points in the figure. Once all elements were computed, we calculated the average of 
the percentage difference between the values of a before and after the median filter was 
applied. If the percentage difference was larger than 5%, we repeated the procedure until 
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convergence was assured. To test that our median filter is actually removing only noise and 
not actual data, we devised the following test. We applied our median filter to one of our 
low-resolution surveys, until the 5% convegence criterion was fulfilled. We then computed 
the difference cTjj between filtered (Fij) and raw data (Rij), and computed the standard 
deviation of the difference a. Wc then generated a fictitious matrix ai/ of noise following 
a Gaussian distribution with standard deviation a and zero mean, and added this fictitious 
noise to the filtered data, obtaining a matrix of a values R[j. The median filter was then 
reapplied to these new data until convergence was assured, and a new matrix F-j obtained. 

To compare the two distributions Fij of filtered data and F-j (data obtained after the 
filtration of the fictitious noise), we devised the following scheme. We computed a variable, 
whose expression is given by: 

where aij and ai/ are the errors on the values of the filtered data, intended as the difference 
between raw and filtered data. The variable so computed is supposed to follow a - 
like distribution. The problem is to determine the number of degrees of freedom of the 
distribution. This is not just given by the number (399) of data points in our low-rcsohition 
survey since, when taking the median, each value of a is connected to the value of at least 
some of its eight neighbors. Moreover, the process of filtering is apphed several times, until 
convergence is assured, so that the number of degrees of freedom can be substantially smaller 
than the number of data points, in a way that is very difficult to estimate. 

However, as a rule of thumb, a value of considerably smaller than 399 is a relatively 
good measure of a good fit. In the case of our test, the value of was 55, so we believe our 
test should be a convincing proof of the validity of our application of the median filter. 
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ro [degrees] m [degrees] 

Fig. 1. — Energy levels of the Kozai Hamiltonian for O = 0.70 and = 0.25, in the uj vs. 1 — 
plane. Since the levels are symmetric about uo = 180°, we only show the interval — 180°. 
For < 0.6 librating solutions are possible (see the right panel). 
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Fig. 2. — The fate of orbits located near the separatrix between circulation and libration for 
our low-resolution survey. Librating particles are shown by full dots, circulating orbits are 
shown by black points; and asterisks denote orbits that switched from libration to circulation. 
The black line identifies the separatrix according to the secular model for S/2000 S5 Kiviuq's 
B at the maximum eccentricity. Kiviuq's position is indicated by a cross (S5 in the figure). 
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Fig. 3. — Contour plots that show the periods of hbration or circulation in uj in years. The 
argument of pericenter lo) lies on the horizontal axis, and on the vertical axis x = 1 — e^. 
The solid parabolic line identifies the separatrix according to the secular model for S/2000 
S5 Kiviuq's B at the maximum eccentricity. The dashed line shows the approximate location 
of the actual separatrix (see Fig. 2). 
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Fig. 4. — Low-resolution survey, grayscale plots of a values for AT =1 Myr. The solid 
parabolic arc in the figure's center represents the separatrix according to the secular model 
(see Fig. 2). Chaotic orbits are associated with lighter shades (i.e., higher values of logio{a), 
(see grayscale at right). 




Fig. 5. — The level of chaos for a simulation where particles are subjected to the gravitational 
perturbations from Saturn and Jupiter. Chaotic orbits are in white, while orbits of regular 
behavior are in black (see grayscale at right). 
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a [degrees] ro [degrees] 



Fig. 6. — Measures of chaos for the simulation with the SEC Brctagnon solution (upper 
left), for the A21SEC Bretagnon solution (upper right), the GISEC solution (lower left), and 
the A21GISEC solution (lower right). Light shades correspond to logio{a) = —3 or to more 
chaotic orbits. 
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m [degrees] m [degrees] 



Fig. 7. — Lyapunov time for integrations (left) when the traditional version of the Lyapunov 
integrator is used and (right) when the SEC solution of the Bretagnon model is applied. The 
left figure has a feature of long Lyapunov times, associated with a region of high frequency 
dispersion, but on the right figure we see that this feature disappears when the SEC solution 
is used b). The substantially different gray scales used in the two figures were chosen to 
maximize the contrast between nearby features of T^. 
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Fig. 8. — Time evolution of the angle zu — — 5Xs + 2Xj + S^ls for an orbit in the region 
of the secondary resonance involving the Great Inequality {xq = 0.756, cuq = 90°). Note how 
the resonant argument alternates between periods of libration and periods of circulation. 
This particle has a large value of a. 
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0.65 0.7 0.75 



x=(1-e') 

Fig. 9. — A plot of logiQ{cr) versus initial x for our high-resolution survey. The dashed vertical 
lines show the transition region in which particles may switch behavior from circulation to 
libration. The two major sources of chaos for the system are this transition region and the 
pericentric secular resonance at x = 0.68. Especially important in the context of planetary 
migration is the secondary resonance between lj and the Great Inequality at x = 0.696. 
Other secondary resonances are also present. 
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t[yr] xio' t|yr] 



Fig. 10. — Resonant arguments for test particles in the a) pericentric secular resonance, b) 
resonance of argument zu — Q — 5\s + 2Xj + 3fls, c) resonance of argument 4:{zu — Q) — 3{5\s — 
2Aj) + 8Qs + cind d) 3{zu — Q — 5Xs + 2Aj + 3Qs) + ~ ^s)- AH these particles shows 
periods of libration for their respective resonant angles. The orbital elements of planets and 
satellites are computed with respect to the solar system's invariable plane. 
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Fig. 11. — Periods of cu's precession obtained with FAM as a function of x. Vertical lines 
represent the locations of features of high cr, the dashed vertical lines show the transition 
region in which particles may switch behavior from circulation to libration, and the horizontal 
lines report the location of the commensurabilities between lo and the Great Inequality. The 
positions of the actual resonances slightly differ from that of the commensurabilities (see 
Eq. 9). 
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Fig. 12. — The period of the Great Inequahty as it varies with positive expansions in Jupiter 
semimajor axis, assuming Saturn's orbit fixed. 
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Fig. 13. — Shaded plots of a when the Great-Inequality's period was (left) 640 yrs, and 
(right) 480 yrs. The position of the secondary resonance shifts upward toward the center of 
the libration island in the second case (480 yrs is the minimum period of libration) . 
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Fig. 14. — Variations in the frequency with largest amphtude as a function of x. When a 
resonance is encountered, frequency values are scattered {x = 0.68 — 0.69). The curve is not 
continuous at x = 0.71 because different frequencies have their largest values of amplitude 
for different regions in x. The dashed vertical lines show the location of the transition region, 
where orbits switch behavior from circulation to libration. 
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Fig. 15. — Variations in the frequency having the largest amphtude as a function of x, for 
three regions in x, and three values of the Great-Inequality's period (810, 883, and 950 
yrs, respectively). The first row of figures reports the region around the 1:1 GI resonance, 
the second the region around the 3:4 and resonance, and the third the region of librating 
particles. Vertical lines show the expected position of the resonance (identified at the top of 
the line) based on its resonant argument. By 3(1:1) we indicate the resonance of argument 
3(w — VL — 5As + 2A J + 3^5) + (f2 — Vts)- Dashed vertical lines in the third row of figures 
show the position of the transition region. See text for more details. 
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Fig. 16. — The time evolution of the semimajor axes of (left) Jupiter and (right) Saturn. The 
black band plots the results of our integrator, the triangles show the result of the integration 
of Eq. 10, and the line displays the results of a simulation on planetary migration from 
Beauge et al. 2002, as described in the test. 
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Fig. 17. — Fate of test particles for our low-resolution survey, where we used Eq. (10) to 
simulate planet migration, a) gives the orbital nature at the simulation's beginning, while b) 
shows it at t = 5 Myr; the symbols are the same used in Fig. 2. At the end of the simulation 
only one additional particle was lost from the libration island. 
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Fig. 18. — Time evolution of uj for an orbit initially inside the Kozai resonance (a;o = 
83,a;o = 85°). At t = 3.7 Myr the particle is captured into the secondary resonance with the 
Great Inequality. As a consequence, the amplitude of libration slowly increases until the test 
particle reaches the separatrix of the Kozai resonance. At that point, it escapes from the 
resonance and starts switching back and forth from circulation to libration (around either 
90° or 270°). At the end of the integration (t = 10 Myr) the particle was on a circulating 
orbit. 
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Table 1: The number of terms in a, X, {K, H), {P,Q) used for four separate solutions of 
Saturn's orbit based on the Bretagnon model. 



Model 


a 


A 


{K,H) 


{P,Q) 


SEC 


1 


12 


2 


2 


A21SEC 


12 


12 


2 


2 


GISEC 


1 


12 


8 


6 


A21GISEC 


12 


12 


8 


6 
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Table 2: Results of a test that compares the fihered results of four Bretagnon simulations 
with the filtered results for a full simulation involving Saturn and Jupiter. Smaller values of 
indicates a better fit (the maximum number of degree of freedom is 399, see appendix). 
The models are described in the text. 



Model 


SEC 


A21SEC 


GISEC 


A21GISEC 




2741.6 


132.3 


70.7 


33.5 
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Table 3: Values of the asymmetry coefficient for four simulations with different initial con- 
ditions, corresponding to particles all having the same values of a) Q — = constant, b) 
- ^0 = 0°, c) - ro© = 180°, A - A© = 180°, and d) 2(ro - A©) = 45°. Case b) seems 
to be the most symmetric, but all results overlap to within the errors. 



Nodal Resonance 
Q — Qq = const. 


Pericentric Resonance 

ZU — ZUq = 0° 


Pericentric Resonance 
w-Wq = 180° 


Evcction Inequality 
2{zu - A©) = 45° 


60±13% 


55±12% 


64±14% 


60±18% 
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Table 4: Resonant arguments for the secondary resonances we observed or suppose to exist 
for the region near the separatrix of the Kozai resonance. The second column reports if 
the resonant argument was observed to librate. We used the word "suspected" or "strongly 
suspected" for resonances inside the libration region, whose positions, when the period of 
the Great Inequality was changed, moved according to our predictions. Other resonances 
involving higher commensurabilities between u {= w — Q) and the Great Inequality are 
weaker, and are not reported. 



Resonant argument 


Observed 


zu-Q-5Xs + 2Xj + 3Qs 


yes 


- - 5A5 + 2Aj + 2zus + 


yes 


w - Q - 5Xs + 2Xj + 2-njj + Qs 


no 


S{w - Q) - 2(5A5 - 2Xj) + bfls + 


no 


3{w -Q)- 2{5Xs - 2Xj) + bwg + 


no 


3(w -9)- 2(5A5 - 2Aj) + bwj + Vis 


no 


A{w -Vl)- 3(5A5 - 2Aj) + 8Q9 + 


yes 


A{w -n)- 3(5A5 - 2Aj) + 9^5 


yes 


A{w -n)- 3(5A5 - 2Aj) + 9u7j 


no 


3{w-Vl- hXs + 2Aj + Sas) + {9- 9s) 


yes 


?>{w - Q - 5A5 + 2Aj + 3^5) + (fi - wj) 


no 


3{w - Q - 5A5 + 2Aj + 3wj) + (Q - vos) 


no 


2{w - Q) - 3(5A5 - 2Aj) + 8Q5 + ws 


strongly suspected 


2{w - fi) - 3(5A5 - 2Aj) + 9vjs 


suspected 


2{w -9)- 3(5A5 - 2Aj) + 9wj 


suspected 


?>{w -9)- 4(5As - 2Aj) + 119s + 


strongly suspected 


3(w -9)- A{hXs - 2Xj) + llws + 9s 


suspected 


3(w -9)- 4(5A5 + 2Aj) + llwj + 9s 


suspected 



